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ABSTRACT 

Crossmatching catalogs at different wavelengths is a difficult problem in astronomy, 
especially when the objects are not point-like. At radio wavelengths an object can 
have several components corresponding, for example, to a core and lobes. Consider¬ 
ing not all radio detections correspond to visible or infrared sources, matching these 
catalogs can be challenging. Traditionally this is done by eye for better quality, which 
does not scale to the large data volumes expected from the next-generation of ra¬ 
dio telescopes. We present a novel automated procedure, using Bayesian hypothesis 
testing, to achieve reliable associations by explicit modelling of a particular class of 
radio-source morphology. The new algorithm not only assesses the likelihood of an 
association between data at two different wavelengths, but also tries to assess whether 
different radio sources are physically associated, are double-lobed radio galaxies, or 
just distinct nearby objects. Application to the SWIRE and ATLAS CDF-S catalogs 
shows that this method performs well without human intervention. 
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1 INTRODUCTION 


To assemble time-series and multicolor data sets 
in astronomy, we routinely combine observations 
from different telescopes and instruments. Hundreds 
of catalogs are available today starting with the 
large st surveys ( e .g., SDSS SkyServer; lYork et ^ 
kOOd : ISzalav et al.l 1200 j ) , to the many legacy ob¬ 
serva tions in Vizier (lOchsenbein. Bauer fc Marcout I 
kOOd l , which are served by Virtual Observatory 
(VO) services as well as Simbad llWenger et al.l 
200(J) or the NASA Extr a galact ic Database (NED; 
Mazzarella. Madore fc Melon I koOll) . NASA/IPAC In¬ 


frared Science Archive (IRSA) Q. These collections are 
accessible online and there are also specialised services 
that offer cross-identifi c ation . Such services include the 
SkyQ uery dMalik et alJ [200^ : iBudavari. Dobos &i Szalav I 
kOld) . CDS XMatch 0 or the VAO Cross-Comparison 


* E-mail:dfan2@pha.jhu.edu 
t E-mail:budavari@jhu.edu 
^ IRSA http://irsa.ipac.caltech.edu/ 

^ CDS xMatch service http://cdsxmatch.u-strasbg.fr/xmatch 


Tool 0. These services use fast indexing meth- 
ods, e.g., the Hierarchical Triangular Mesh (HTM; 


^^— ff 11 ° 

Kunszt. Szalav & Thakai 20011: Budavari. Szalav 

V Fekete 1 

2010l). HEALPix (Gorski. Hivon & Bandav 


20051) 

or the Zones Algorithm (iGrav et al. 

2004|; 


ICrav, Nieto-Santisteban fc Szalavl l200d ) that quickly 
eliminate sources too far away from each other. 


These tools, however, are essentially designed for cross¬ 
matching between point sources. They are good for merg¬ 
ing optical catalogs of stars and small extended sources but 
mostly useless for associating the largest galaxies or radio 
sources. For example, a single radio galaxy can have sev¬ 
eral distinct components corresponding to its core and lobes. 
Traditional methods can quickly link the radio core to the 
optical center, but the jet or lobes would be ignored since 
they are far away from the center. Currently such compo¬ 
nents are identified and associated by eye (e.g.. lNorris et ahl 

l2006l) . 

Matching by eye is adequate for existing deep radio sur¬ 
veys, with only a few thousand sources, but will be untenable 


^ VAO Catalog Cross-Comparison Tool http://vao- 
web.ipac.caltech.edu/applications/VAOSCC/ 
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for some surveys currently under way with tens of thousands 
of source s, and will be tota lly impractical for surveys such 
as EMU (|Norris et al.ll2oill ~l which is expected to detect 70 
million sources. Existing attempts to build automated al¬ 
gorithms to perform the cross-matching of complex radio 
sourc es fe.g.. Ilhoctor 1120061 : Ivan Velzen. Falcke fc Kordingl 
l2015h . have used simple algorithms that will only identify 
the simplest cases. 

Here we present an automatic method to estimate the 
matching probability for radio sources, which e xtends the 
Bayesian approach of iBudavari fc Szalav I (l2008ll to include 
an explicit geometric model. Eor the current analysis, this 
model is aimed at simulating a typical radio galaxy struc¬ 
ture, composed of a core and two approximately symmet¬ 
rically located lobes. In principle, though, any model can 
be used, and future work will generalise this to account for 
more complex radio source structures. Our approach is de¬ 
scribed in section [21 Its application and results are discussed 
in section [3] Section [4] presents our conclusions. 


2 MATCHING WITH RADIO MORPHOLOGY 

Probabilistic cross-identification has several advantages over 
traditional methods. Bayesian hypothesis testing uses in¬ 
dividual estimates of the uncertainties in the directions of 
sources and provides an objective measure of the quality of 
a given association from the available data. It also just re¬ 
quires a single threshold value to define a matched catalog 
even if there are many input catalogs to consider. Its biggest 
advan tage, however, i s its fl exibility to model the problem at 
hand. iKerekes et al.l (I 2 OI 0 I I successfully applied the method 
to associate mov ing stars with unknown proper motion and 
iBudavain (|201ll ) extended the method to time-domain ob¬ 
servations to find matching events, e.g., supernovae. Here 
we use the same concept to associate radio sources with an 
unknown number of components and geometry. 

To decide whether a given association is a good match, 
the likelihoods of two hypotheses have to be calculated. In 
practice these are done as integrals over the parameter space 
of the model. For point sources this is a simple ca lculation 
with an analytic result (iBudavari fc Szalav Il200a 'l. but in 
general it would have to be computed numerically. 

Throughout this paper we use the flat-sky approxima¬ 
tion local to each association. The point source’s nearby 
coordinates are projected to its tangent plane. This opera¬ 
tion converts the spherical geometrical problem to a classical 
plane geometrical issue. This approximation is very accurate 
for the small distances we consider and allows for fast cal¬ 
culations. 


2.1 Choosing a prior 

For this analysis we assume that a radio source complex is 
well modelled either by (in the case of a star forming galaxy 
or radio-quiet AGN) a single isolated object , or (in the case 
of a radio-loud AGN) by a radio nucleus, coincident with 
the host galaxy, together with a pair of lobes , resulting in 
a linear or near-linear arrangement of 2 or 3 components. 
In some AGNs, either the nucleus or one of the lobes may 
be missing or below our detection limit. In the absence of a 
detailed quantified distribution of radio source morphologies 


this is the simplest prior to assume, and appears to be valid 
for the vast majority of known radio sources. 

A more comprehensive prior could be constructed from 
the distribution of radio sources found in existing sur - 
vey d ata, such as FIRST! Becker. Wlute_^£jIelfaJid_ 


NVSSlICondon et al.lll998ll . SUM SSlfeock. Larp fc Sadler I 


I 1999 I : iMauch et al. 2003ll . ATLAS dNorris et al.l 2006h . This 
distribution would vary as a function of sensitivity and angu¬ 
lar resolution, and so the construction of such a prior would 
be justihed only if the simpler approach taken here proved 
inadequate. However, we recognise that such a comprehen¬ 
sive prior may be required in subsequent implementations 
of this technique. 


2.2 Straight-Line Model 

As an example, we consider a simple triple radio source 
model, consisting of a radio core and two ejected lobes, 
equidistant from the core, with all three components lying 
along a straight line. The measured positions of the core 
and two lobes are denoted by the vectors m, m' and m" 
respectively, as shown in Fig. [T] These are two-dimensional 
vectors on the tangent plane and, in this simplest case, obey 
the following relation 


m = 


m' -I- m” 
2 


( 1 ) 


which yields m” = 2m — m' as illustrated in Figured] 

We must now compare this model with real data, con¬ 
sisting of the measured position of a putative optical host 
galaxy denoted by c, and the measured positions of an as¬ 
sumed radio core and two lobes, denoted by ro, ri and r 2 
respectively. We label the measurements such that c and 
ro correspond to the core at m, while ri and r 2 are the 
measurements of the lobes at m' and m”. If the detections 
come from the same object with the core and lobes iden¬ 
tified correctly, the likelihood of this hypothesis is calcu¬ 
lated as the integral over m and m'. The prior on m is just 
p(m) = l/47r as the object can be anywhere on the entire 
sky but the conditional probability density of m' is local¬ 
ized around m, e.g., p{m'\m) = /(|m' —m|). We will study 
the effect of this choice later. In this notation, the likelihood 
of the CORE-LOBE-LOBE configuration is 


C-cll ~ 


Jdm p{m) Lc{m) Lrgim) 

Jdm' p{m'\m) Ln (m') Lrj {2m —m') 


( 2 ) 


where the likelihood functions have the usual form, say a 
Gaussian Lx{m) = G{x\ m, a^) with the appropriate uncer¬ 
tainties for each detection. This is a two dimensional dou¬ 
ble integral including all the components in the geometrical 
model. Lobes depend on the core so the core’s position m 
needs to be established first, followed by the lobes m! and 
m". 

Naturally, there are many other competing hypotheses, 
where we consider only partial matching. For example, we 
consider a LOBE-LOBE model in which the putative core 
is actually an unrelated radio source, we write the LOBE- 
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LOBE likelihood as 


Cii = 


Idmo p{mo) Lrg (mo) 


Jdmi p{mi) Lcim-i) 


jdm'i p(mi |mi) (mj) (2mi — m[) 


(3) 


where the first term can be evaluated analytically in case of 
the fiat all-sky prior on mo and is equal to l/dvr because 
the Gaussian is symmetric in its arguments. The first term 
is independent of the second, as their host galaxies are not 
the same, s o their integrals are calcu lated separately, then 
multiplied (iBudavari fc Szalav Il2008h . 

In practice there can be more than three radio sources in 
close proximity to a given detection in another band, which 
are considered to be separate sources like the mo. General¬ 
ising our model to such complex scenarios will be explored 
in future work. 

A more realistic model would also allow the lobes to be 
at slightly different distances. We introduce a new parameter 
k to account for such asymmetry and write 

m" = 2m — m' + k{m — m'). (4) 

The numerical integration now includes the k parameter 
over a prior that we choose to be a Gaussian centered on 
zero, so fc = 0 corresponds to the symmetric case. 

Each of the radio components can be considered a core, 
a lobe or a separate object entirely. For example, three 
radio components (ro,ri,r 2 ) can be CORE-LOBE-LOBE, 
LOBE-CORE-LOBE, LOBE-NONE-LOBE, NONE-LOBE- 
LOBE, LOBE-NONE-CORE and so on. We calculate the 
likelihoods of all these possibilities. The ratio of possible 
likelihoods for two hypotheses with likelihoods C\ and C 2 

B 1.2 = (5) 

is the Bayes Factor that compares them in an objective 
way. When BS>1, the data clearly prefer the hypothesis 
in the numerator but when Bi ,2 ~ 0 the alternative has 
stronger support. When Bi ^2 « 1 is around unity, more ev¬ 
idence might be needed to make a decision. 


2.3 Implementation Details 

Our catalogs are stored in a database indexed by the 
Zone Algorithm (iGrav. Nieto-Santisteban fc Szala^ l2006l i 
to quickly identify possible candidate matches and rule out 
the completely unlikely ones. We consider lobes with angular 
separations up to 2 arc minutes. These associations are pro¬ 
cessed by a custom multi-threaded C# library that imple¬ 
ments the combinatorics and the numerical integration. For 
a given set of detections we generate all possible hypothe¬ 
ses using the Facet. Combinatorics package by AkisonQ and 
prune the resulting list to only contain the physically pos¬ 
sible combinations, e.g., no more than one core. We further 
exclude models that will clearly perform poorly based on 
the vicinity of the radio sources to the point source in the 
other band. This filtering, however, is very conservative and 
the hnal decision is based on the probabilistic analysis. 


To scale to the large number of possibilities, the numeri¬ 
cal integration has to be very fast. Our impl ementation uses 
impo rtance sampling for all integrals (e.g., I Flannery et al.l 
Il992li . We speed up the calculations by always sampling 
from the tightest possible Gaussian. The product of Gaus- 
sians and linear parameter transformations yield Gaussians, 
which makes the procedure an order of magnitude faster 
than naive implementations. 

The models can contain at most three radio compo¬ 
nents to fit a core and two lobes. When the number of 
components is more than 3, e.g., there are 5 nearby radio 
detections, we consider the extra ones to be from other in¬ 
dependent sources. For example, the NONE-LOBE-LOBE- 
NONE-GORE hypothesis, (where we use the “NONE” des¬ 
ignation to indicate that the respective radio component is 
not related to the counterpart of interest), will have a likeli¬ 
hood that is the product of an integral identical to the previ¬ 
ous Cell formula with the appropriate data and two indepen¬ 
dent integrals that are again just I/Att. For the all-NONE 
hypothesis corresponding to 5-1-1 independent sources, the 
likelihood is a constant (l/47r) . We can compare all other 
hypotheses to this, including those with a different num¬ 
ber of components, by considering additional independent 
sources. 


3 APPLICATION TO SWIRE AND ATLAS 

We apply the new algorithm to two well-studied catalogs. 
The Australia Telescope Large Area Survey (ATLAS; Norris 
et al. 2006, hereafter Norris06) detected 784 radio compo¬ 
nents from 726 distinct sources around Ghandra Deep Field- 
South are a. The Spitzer Wide -Area Infrared Extragalactic 
fSWIRE: iLonsdale et al.]|2003li contains about 113000 ob¬ 
jects in the area also covered by ATLAS GDF-S. 

We compare our associations to the matched catalog 
bv iNorris et al.l ll2006li . which is publicly available online |j 
These associations were created by expert examination of 
the radio contours overlayed on the images, as illustrated in 
Figure [2] There are 119 SWIRE objects that are not avail¬ 
able publicly, so we have removed them from further compar¬ 
ison. The final public sample contains 10 radio triplets and 
27 do ublets. Some doublets are identified by I Norris et al.l 
(l2006li as two lobes, while others are identified as a lobe 
and a core. There are also 558 single sources, which corre¬ 
spond to a single SWIRE source in their list of associations. 
Naturally we do not know the underlying physical reality of 
the published associations but use these published results as 
our reference, against which we benchmark our automated 
method. 


3.1 Most Likely Associations 

To find associations, first we look for the potential radio 
components around each SWIRE detection using an inclu¬ 
sive search radius of 2 arc minutes. This heuristic limit is 
motivated by previously identified associations in the refer¬ 
ence catalog. Its actual value does not affect the quality of 


Permutations, Combinations, and Variations using C# Gener¬ 
ics http://www.codeproject.eom/Articles/26050/Permutations- ® SWIRE & ATLAS CDF-S matched catalog http://cdsarc.u- 

Combinations-and-Variations-using-C-G strasbg.fr/viz-bin/Cat?cat=J/AJ/132/2409 
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Table 1. Comparison between the result of automated and man¬ 
ual classification 


Type 

NorrisOG 

Bayesian 

Common 

Missed 

Extra 

triplet 

10 

12 

9 

1 

3 

doublet 

27 

48 

19 

8 

29 

single 

558 

566 

536 

22 

30 


the automated associations as long as it finds all candidates, 
but a large threshold would slow the algorithm down simply 
because of the large number of candidates to sort through 
and reject. Next we evaluate the likelihoods for all possible 
combinations. In principle, this is combinatorially expensive 
and would require a careful analysis of prior probabilities 
of the different types of matches. In this study, we apply a 
simpler algorithm instead to look for the most likely associ¬ 
ations. 

We use a simple approach to pick the associations based 
on the calculated likelihoods for the different hypotheses. 
It is an iterative procedure where we take the association 
with the largest likelihood. Next we remove from the fist 
of available hypotheses all of those that include the source, 
which we already accounted for in the best association. The 
remaining hypotheses are again searched for the best associ¬ 
ation and we repeat the procedure until nothing is left. This 
greedy approach runs quickly and provides robust results in 
our tests. 

To compare the different hypotheses one has to char¬ 
acterize the possible geometric arrangements of the cores 
and lobes. In particular the probability distribution function 
(pdf) of possible core-lobe distances is needed for the prior. 
This appears as the conditional density p(m'\m). This func¬ 
tion is expected to rise at small separations but fall at large 
scales. One possible such pdf is the Rayleigh distribution. 
In particular we pick a pdf with a mean of 9 arc seconds, 
chosen as a reasonable value that is consistent with the dis¬ 
tribution of realistic core-lobe separations in radio images 
of the sensitivity of ATLAS. Table [T] illustrates the quality 
of our results. At first glance we see that the automated 
method using only the coordinates and a simple geometric 
model can find most of the triplets in our references catalog 
and successfully identify the singletons, too. The doublets 
are associated with less precision: while the new method 
misses some of the reference matches, it also identifies extra 
ones based on the limited information included. In the next 
section we will discuss these categories in detail and explore 
examples to illustrate the scenarios where the reference and 
automated catalogs differ. 

3.2 Association of Radio Triples 

The automated code finds 9 of 10 of the manually identi¬ 
fied triplet associations. Figure [3] shows some of the triplets 
that were common across the reference and automated cat¬ 
alogs. We see that there are a variety of shapes that are not 
only often asymmetric in the distances of the lobes from the 
core but also have a slight angle from the assumed straight 
model. Our method finds these because the astrometric un¬ 
certainty allows sufficient deviation from the straight line. 


These results have high likelihoods on an absolute scale and 
also compared to the possible competing scenarios. 

The only missed triplet has significantly lower likelihood 
compared with the others. Figure [4] shows this association. 
The difference is obvious upon inspection: this geometry is 
more bent and has an angle of about 15° from the straight 
line. This appears to be too much to be tolerated by the 
simple model. A possible improvement of our geometric de¬ 
scription would be to include another parameter to capture 
the bent shape. We plan to explore this possibility in future 
studies. 

Our program also finds 3 extra triplets, seen in Figure[Sl 
Here we discuss these in detail, to guide future refine¬ 
ment of the algorithms. 

C761 

The algorithm classifies C760, C759, and C761 as a ra¬ 
dio triplet. NorrisOG classified them as a “linear complex 
(jets or arc?)”. Both NorrisOG and the algorithm associate 
this source with the SWIRE source SWIRE3-J033533.90 — 
273310.9, and so the only difference is one of terminology, 
caused by the complex shape of the radio source. Classifying 
such complex structures is beyond the scope of the present 
project. 

C349 

The algorithm classifies C351, C349, and C345 as a radio 
triplet. NorrisOG classified them as a “a line of three galax¬ 
ies”. C351 and C349 are each associated with a catalogued 
SWIRE galaxy, whilst C351 is associated with an uncata¬ 
logued SWIRE source. Given that each of the three radio 
components is associated with a SWIRE source, NorrisOG 
classified them as there separate galaxies, whereas the al¬ 
gorithm, without access to the uncatalogued data, found 
matches for only two of the galaxies, and so classified them 
as a triplet. If the algorithm had had access to the same 
catalog as NorrisOG, then it would have made the same clas¬ 
sification. 

C371 

The algorithm classifies CG70, C671, and CG74 as a radio 
triplet. NorrisOG classified them, along with nearby CG75, as 
“A line of 4 galaxies”. CG70 and CG74 are each associated 
with a catalogued SWIRE galaxy, whilst CG71 and CG75 are 
associated with uncatalogued SWIRE sources. Again, if the 
algorithm had had access to the same catalog as NorrisOG, 
then it would have made the same classification. 

The lesson from this is that the robustness of the cross¬ 
matching algorithm, is limited by the availability of mul¬ 
tiwavelength data. NorrisOG were able to use the SWIRE 
images which show sources that were absent from the cata¬ 
log. Modifying the algorithm to enable access to such lower 
reliability multiwavelength data is deferred to future work. 

3.3 Association of Radio Doubles 

Among the 19 doubles that are common in the reference 
and automated catalogs, 14 are LOBE-LOBE combinations. 
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Three typical examples are shown in Figure [6] There are 
eight doubles that our code does not find. In four of these 
cases the automated code associates the radio pairs with 
other SWIRE objects, in combinations that appear to have 
larger likelihoods. In the four remaining cases, the code splits 
the radio components, and rearranges them with other radio 
detections. Figure [3 shows one of these misidentified associ¬ 
ations, which we now discuss in detail as a case study. 

The algorithm associates C473 and C478 as a double, 
and classifies C477 and C476 as isolated galaxies. NorrisOG 
identify C477 and C478 as a pair of galaxies, and C473 and 
C476 as a double-lobed radio source, because there is a low¬ 
brightness bridge of radio emission between C473 and C476, 
so these are likely to be associated, and between C477 and 
C478, so these are likely to be associated. Furthermore, C477 
and C478 each have a bright IR counterpart, so are likely 
to be galaxies. Neither C473 nor C476 have an IR source 
coincident with the radio peak, but there is a bright IR 
galaxy part-way along the bridge, which is likely to be the 
host, especially since the bridge is bright at that position. 

To improve the algorithm, greater weight needs to be 
given to the presence of IR counterparts (such a C477 and 
C478). More importantly, visual classification is advantaged 
by the potential to make use of extended or low signal-to- 
noise features, such as the faint extensions from the radio 
source peaks that link the lobes in these examples. This in¬ 
formation is not captured by the simple positional informa¬ 
tion retained in a source catalog, limiting the performance 
of the automated algorithm in this situation. It is difficult to 
see how such information can readily be incorporated into a 
Bayesian algorithm. 

Among the 29 extra doubles found by the automated 
procedure, there are 14 LOBE-LOBE combination and the 
rest are all CORE-LOBE. While many of these may be 
real, some are clearly two radio components sitting within 
an extended SWIRE object, resulting their LOBE-LOBE 
or CORE-LOBE hypotheses having higher likelihoods than 
being isolated sources. Although this leads to the classifica¬ 
tion of such sources as radio doubles, they are unlikely to be 
so in reality. Ensuring the automated process can recognize 
and exclude this kind of false-double is part of the planned 
refinements. 


4 CONCLUSION 

We have introduced a novel approach to automatically cross¬ 
match radio sources having a realistic morphology against 
point source counterpart catalogs. Instead of the usual cata¬ 
log matching, where detections in separate datasets are sim¬ 
ply assigned to each other when co-located, we can also iden¬ 
tify components of the radio sources as cores or lobes sepa¬ 
rated from each other on the sky. Traditionally such match¬ 
ing could only be done manually by astronomers through 
visual inspection of the overlays of images. 

The new method was tested on the SWIRE infrared 
and ATLAS radio CDF-S catalogs. We find that the au¬ 
tomatic associations agree well with the reference catalog. 
The differences all arise either from inadequate representa¬ 
tion of the data by the catalog, or from the simplicity of 
our current radio source morphology model, which only al¬ 
lows for the components to lie along a straight line. Also the 



Figure 1. A simple model of the radio morphology has a core 
and two lobes. The model coordinates are compared with the 
detections using the usual likelihood function based on a Gaussian 
astrometric uncertainty. 

current procedure only uses the directions of the measured 
sources, while experts constructing a manual catalog would 
typically rely on shape and brightness information as well. 
This study is a first step toward automating a manually 
intensive process in order to provide a scalable functional¬ 
ity for the coming generation of massive radio surveys with 
Square Kilometre Array pathfinder telescopes. 
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Figure 3. An array of sample triplets found both by Norris et al.(2006) and by our method. Each of these triplets has a central core 
and two lobes nearly in a straight line. 
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Figure 5. Three extra triplets found by the program. Most of them are fit well with the geometry model, but are identified as a line of 
galaxies or other cases by Norris et al.(2006), see text for detailed discussions. 
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Figure 6. Three typical radio doubles found by both methods. The left figure is a CORE-LOBE combination. The middle and the right 
figures are of LOBE-LOBE type. 



Figure 2. A typical cross-identification from Norris et al.(2006). 
Contours show the radio emission, and greyscale show the SWIRE 
3.6 fim emission. The three radio components (C634, C631, C628) 
were identified by Norris et al.(2006) as a radio triple source 
hosted by the bright infrared galaxy coincident with C631. 
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Figure 4. A triplet with a large angle between the two lobes. Our 
program classifies it as a LOBE-LOBE which has larger Bayes 
Factor comparing to the triplet. 



Figure 7. A radio double that our program did not find. The al¬ 
gorithm identified the radio objects C478 & C473 as radio double, 
and this association has a larger Bayes Factor than the correct 
pairs, as discussed in Section 13.31 
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